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. In the thermodynamic hmit, a probabihstic cellular automaton can be approximated by 

' a deterministic nonlinear map. Here we construct such a map for the forest fire problem. 

(N ■ The construction is based on the results of the Monte Carlo simulation, performed on 

a square lattice of million cells. The results of the calculation are analyzed by means of 
the Hoshen-Kopelman algorithm (HKA). The only parameter of the map describes the 
probability that a tree appears at an empty cell during one time step. The obtained map 
seems to be non-differentiable at the percolation threshold. The Lyapunov exponent for 
' the map is positive. Also, we found the cycle of length three by means of the method 

, of symbolic dynamics. The results are illustrated by the experimental data on the 

' forest fires in Canada in years 1970-2000. Although these data are fortunately far from 

I thermodynamic limit, their qualitative character is reproduced for smaller lattices. 
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^ ! 1 Introduction 

\ The problem of forest fires has got a wide audience [0 for its percolation origin |^ and rich 

O ■ connections to basic problems in computational physics, e.g. to the self-organized criticality. 

On the other hand, the problem is not entirely virtual, although experiments are not eagerly 
accepted by the society. The data can be found in annual statistics of forest rangers and 
insurance companies. Both institutions are obviously interested in predicting of forest fires, 
at least in the statistical sense. Hovsfever, there are some indications that such predictions 
I are not possible in long time scale — the problem is chaotic |§, |^ . 

The present work is a continuation of Ref. In that paper, a map p„ — > Pn+i which 
describes the forest density p in subsequent time steps n was proposed. The construction 
was based on the idea of reducing a probabilistic cellular automaton to the map, which is 
equivalent to the mean field approximation of the automaton rules . The map was designed 
not so as to reproduce the dynamics of the fire front, but rather to simulate the statistics of 
the forest density in subsequent years n. The forest was represented by a square lattice Lx L, 
occupied by trees with probability p. The above map was equivalent to a superposition of two 
functions, where each function transforms the forest concentration p. The growing process 
was represented by a parabolic function, Pn+i = Pn + fPnO- — Pn), where r is a parameter. 
Next, a season of fires was modeled by a reduction of Pn+i by what was burned. This amount 
was demonstrated to be equal to the mean squared size of clusters of trees. By definition, 
a tree belonged to a cluster if it was a nearest neighbor of another tree which belongs to it. 
The neighborhood was of von Neumann kind. 

The prescription of calculating the burned area contained the probability that a cluster 
is ignited. This probability can be defined "per area" or "per cluster". We have obtained 
two kinds of statistics, which can be interpreted as follows: In the former case, a lightning 
strikes at a cell, occupied by a tree or not. The probability that the ignited tree belongs to 
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a cluster is proportional to the size of this cluster, i.e. the number of trees in it. The burned 
area is also proportional to this size. Then, the totally burned area is the mean squared size 
of the cluster, as stated above. The weight of the cluster of size Si is = Si/L'^. In the 
latter case, the fire is set intentionally, what means that a tree is ignited with probability 
one. The rest of argumentation is the same as in the former case. The appropriate weight in 
a such situation is 

w': = s,iY,s,. (1) 

i 

The former case {w') was discussed thoroughly in Ref. In the present work we focus on 
the latter case, where the weights are given by Eq. (0). 

The paper is organized as follows: The next section is devoted to the calculation of the 
map Pn — ^ Pn+ii which is obtained by an application of the HKA |Q. In Section 3 we 
demonstrate the results on the Lyapunov exponent and the symbolic dynamics, which allows 
to obtain the cycle of the length three in the time evolution of the trees density p. There, we 
also refer to the experimental data. Our conclusions arc given in the last section. 



2 The model 

The map p„ Pn+i is a superposition of two functions. The first one, p ^ p + rp(l — p), 
describes the season of growing. As a result, trees appear in empty cells with the probability 
rp(l — p), where r is a parameter from the range [0, 1]. Afterwards, the burned area of the 
forest is calculated with the weights w'l ^ i.e. 

i 

and it is subtracted from trees concentration after growth season, i.e. p ^ p — A{p). 

The cluster size distribution is obtained by means of HKA. With HKA § we arc able 
to label all occupied sites on a L x L large square lattice in such a way that the sites with 
the same label belong to the same cluster and different labels are associated to different 
clusters. As a result, we obtain the cluster size distribution for 1001 values of p, pi = 10~'^i, 
z = 0, 1,2, 10'^. The simulation is carried out on a square lattice 10"^ x 10'^ large, with a 
fraction of p + rp(l — p) occupied sites for different values of the initial trees concentration 
< p < 1. The map is got for 1001 equidistant values of p as an average over hundred trial 
lattice configurations for each value of pi. For other values of p we use straight lines joining 
neighboring points pi and Pi+i. The obtained map is shown in Fig. |l| for three different 
values of the parameter r. 

The method above suffers from the finite size of the lattice and the finite grid. In par- 
ticular, we observe some shift of the maximum of the map for ?■ = with respect to the 
percolation threshold for the square lattice. Instead to occur at pc = 0.59273 we find 
the maximum near p = 0.583. The rounding effect of the finite size reduces the map to 
the well-known unimodal map, where the period-doubling bifurcation must occur, at least 
at some range of the parameter r. That is why we extrapolate the shape of the curve to a 
tentative limit of infinite size (the so-called thermodynamic limit). The result for r = is 

\Pn <^ Pn < Pc, 

Pc exp(- 12. 735 VPri -Pc) Pn > Pc- 



3 Results 

The symbolic dynamics is known to be the most rigorous method to investigate chaos in 
one-dimensional maps js). The basic tool is to assign symbolic words to the periodic cycles. 
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w" (provoked), L=1000, ^=100 




r=0.0 
r=0.5 
r=1.0 

p(n+1)=p(n)+rp(n)(1-p(n)) 



Figure 1: The map for provoked fires {w") and different values of the parameter r. 



Example giving, the cycle of length two is represented by CR, of length three — by CRL. 
etc. The appearance of the cycle three in the time evolution of a given variable is known 
to be a proof of the presence of cycles of all possible lengths. In our case the cycles of the 
variable p are not superstable, because the curve given by Eq. (^) is not differentiable at 
p = Pc- In fact, the cycles are even not stable: we do not observe any windows of stability 
in the bifurcation diagram. Still, it is not difficult to find a cycle three by the method of 
"shooting" . Fixing the value of pi , we change the parameter r until p^ = pi- An example is 
the cycle p = 0.340000, 0.419593, 0.505974 for r ^ 0.354694. 

In Fig. ||we show the values of the Lyapunov exponent as dependent on the parameter r. 
These results are obtained for the extrapolated map given by the Eq. and by the method 
described in Ref. |9| (Fig. ^(a) ). To obtain the Lyapunov exponent we have investigated 
the subsequent differences between two initially nearby trajectories. The slope of the time 
dependence of the logarithm of this difference (see Fig. P) gives the exponent A. The results 



for the finite lattice are qualitatively the same (Fig. 2(b)[ 

In Fig. I we show two plots: calculated time dependence of A{p) for the lattice of 50 x 50 



trees (Fig. f(b)D and the statistical data on the forest fires in Yukon |^, (Fig. 4(a)) 



We have also checked that the phase portraits A{pn+i) vs. A{pn) are more or less the same 
as the appropriate phase portraits for the experimental data [|l0| (Fig. ||). As we can see, 
small fires occur more often than large ones. Both kinds of the phase portraits show some 
anticorrelations between the forest fires in subsequent years: large burned area in a given 
year is likely to be followed by a small area in the subsequent year. This remains true also 
for the data on the forest fires in other provinces of Canada in years 1970-2000. 

The bifurcation diagram for the map obtained by HKA (shown in Fig. |l|) presents itself 
as a homogeneous spot (see Fig. ^). The same is true for the bifurcation diagram for the 
analytic map defined by Eq. (||). The histogram of the forest fires for the map obtained by 
HKA is presented in Fig. |7[ together with the histogram for the analytic map. The difference 
is that in the latter case, small fires do not occur. Still however, the events without fires are 
the most frequent in all the cases. Both histograms show large amount of events without 
fires, and the exponential shape in the range of intermediate and large fires. In this region 
N{A) increases with A, on the contrary to the predictions of the models which display the 
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Figure 2: Linear fit to the first fifty points of Fig. fallows to evaluate the Lyapunov exponent 
A for different r and (a) infinitely large or (b) finite lattice. Only positive values of A are 
presented. 
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self-organized criticality This can be interpreted as a consequence of our assumption, 
that large clusters are ignited more likely. The obtained curve is shifted to right when the 
parameter r increases, as it is shown in Fig. 1^ for the analytic map. 



4 Conclusions 

The results presented allow to state that it is not possible to predict the statistics of the 
forest fires in long time scale. This conclusion is drawn from our much simplified model of 
the fires, where the forest is reduced to a lattice of cells, and a tree ignition — to the nearest- 
neighbor relation with a tree which already burns. Obviously, almost all technical details of 
a real fire Q are omitted in this approach. Still we believe that our conclusion is justified. 
The argument is as follows: suppose that a much simplified model predicts chaos. Suppose 
that the model has been improved by introducing several complex details. Is it possible that 
the new version gives a regular and predictable behavior? The answer is: no. The chaotic 
character could be eliminated by a noise, but the predictability cannot be improved fl^ . 
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Figure 4: The comparison of the forest fire statistics: Time series of fire sizes (a) for Yukon, 
Canada, 1971-1998 and (b) the computer simulation on small lattice. 
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Figure 5: The comparison of the forest fire statistics: A{pn+i) vs. A{pn) (a) for Yukon, 
Canada, 1971-1998 and (b) the computer simulation on small lattice. 
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Figure 6: The bifurcation diagram for the map obtained with HKA. 
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Figure 7: The histogram of burned area obtained for the map obtained by HKA (pluses) and 
by means of computer simulation for the map given by Eq. (^) (dots). Subsequent lines of 
dots correspond to larger and larger value of parameter r = 0.1, 0.2, • • • , 0.9 from the picture 
top to bottom. 
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